#!/bin/csh
# timeseries.csh
####################################################################
# Plots timeseries of different fields
#
# - Set the simulation directory containing YACC's output 
#   in the file "simulation.dir"
# - Run the script using:
#   ./timeseries time0 time1
#   where time0 and time1 are initial and final times of the series
#   If they are not given, the whole time series is plotted
#   By default, the timeseries of all fields marked with 1 are plotted
########################################################################

#-----------------------------------------------------------------------------------
# Possible time-series
#-----------------------------------------------------------------------------------
# temperature (default: 1)
# velocity (default: 1)
# nusselt (default: 1)
# viscosity (default 0)
# viscdiss: viscous dissipation (default 0)
# composition (default: 0)
# heats: Rayleigh number for internal heating (default: 0)
# mixfrac: mixed fraction (default: 0)
# depletion (default: 0)
# staglid: stagnant lid thickness (default: 0)

# Select timeseries to plot
set temperature = 1
set velocity = 1
set nusselt = 1
set viscosity = 0
set viscdiss = 0
set composition = 0
set heats = 0
set mixfrac = 0
set depletion = 0
set staglid = 0
#-----------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------
# Options:
#-----------------------------------------------------------------------------------
# dirout: directory of the YACC's output files (needs to be written in the 
#         'simulation.dir' file)
# dirplots: directory where to store the plots (default: $dirout/Plots/Timeseries)
# makejpg: make also jpg files (default: 1)
# jpgdensity: resolution of the converted jpg files (default: 150)
# deletepdf: delete pdf files (done only if makejpg=1, default: 1)
# dimx: plot width in cm (default: 20 cm)
# dimz: plot height in cm (default: 10 cm)
# X0: x-coordinate in cm of the bottom-left corner of the plot (default: 3 cm)
# Y0: Y-coordinate in cm of the bottom-left corner of the plot (default: 5 cm)
# zmin: z-coordinate of the bottom edge of the plot (default: 0)
# zmax: z-coordinate of the top edge of the plot (default: 1)
#
# Select options
set makejpg = 1
set jpgdensity = 150
set deletepdf = 1
set dimx = 20
set dimz = 10
set X0 = 3
set Y0 = 5
#-----------------------------------------------------------------------------------

source simulation.dir

if (! -d $dirout) then
  echo $dirout 'does not exist. Check the path stored in the variable "dirout" '
  exit
endif

set dirplots0 = {$dirout}Plots/
set dirplots = {$dirout}Plots/Timeseries/

#-------------------------------
# Select YACC's output files
#-------------------------------
set timeseries = ${dirout}timeseries.dat

if (! -d $dirplots0) then
  mkdir $dirplots0
endif
if (! -d $dirplots) then
  mkdir $dirplots
endif
echo 'Plot time series for model:' $dirout

#-----------------------------------
# Set temporary ps output file
#------------------------------------
set output = timeseries.ps

#----------------------------------
# Set gmtdefaults
# ---------------------------------
source setgmtdefaults

#################
# Time-series 
#################
set time0 = ${1}
set time1 = ${2}
if (${1} == "" && ${2} == "") then
  set time0 = `awk 'NR > 1 {if(min==""){min=$1}; if($1< min) {min=$1}} END {printf "%12.5f", min}' $timeseries`
  set time1 = `awk 'NR > 1 {if(min==""){min=max=$1}; if($1>max) {max=$1}} END {printf "%12.5f", max}' $timeseries`
endif
set atime = `echo "scale=3; ($time1 - $time0)/5" | bc -l`

####################
# Temperature
####################
if ($temperature == 1) then
  set adaptive_scale = 0
  echo "Temperature time series"

  set title = ""
  set xlabel = "Time"
  set zlabel = "Temperature"

  set pdfoutput = temperature_ts.pdf
  set jpgoutput = temperature_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk 'NR > 1 {if(min==""){min=$3}; if($3< min) {min=$3}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk 'NR > 1 {if(min==""){min=max=$3}; if($3>max) {max=$3}} END {printf "%12.5f", max}' $timeseries`
    if ($min == $max) then
      set min = 0
      set max = 1
    endif
  else 
     set min = 0.2
     set max = 0.8
  endif
  set delta = `echo "scale=4; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba{$atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk 'NR>1 {print $1, $3}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif


####################
# Velocity
####################
if ($velocity == 1) then
  set adaptive_scale = 0
  echo "Velocity time series"

  set title = ""
  set xlabel = "Time"
  set zlabel = "Velocity"

  set pdfoutput = velocity_ts.pdf
  set jpgoutput = velocity_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk 'NR > 1 {if(min==""){min=$4}; if($4< min) {min=$4}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk 'NR > 1 {if(min==""){min=max=$4}; if($4>max) {max=$4}} END {printf "%12.5f", max}' $timeseries`
    set delta = `echo "scale=0; ($max - $min)/5" | bc -l`
  else
    set min = 0
    set max = 8000
    set delta = `echo "scale=0; ($max - $min)/5" | bc -l`
  endif

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba${atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk 'NR>1 {print $1, $4}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

####################
# Viscosity
####################
if ($viscosity == 1) then
  set adaptive_scale = 1
  echo "Viscosity time series"

  set title = ""
  set xlabel = "Time"
  set zlabel = "log@-10@-(Viscosity)"

  set pdfoutput = viscosity_ts.pdf
  set jpgoutput = viscosity_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk 'NR > 1 {if(min==""){min=log($6)/log(10)}; if(log($6)/log(10)< min) {min=log($6)/log(10)}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk 'NR > 1 {if(min==""){min=max=log($6)/log(10)}; if(log($6)/log(10)>max) {max=log($6)/log(10)}} END {printf "%12.5f", max}' $timeseries`
    if ($min == $max) then
      set min = -1
      set max = 1
    endif
  else
    set min = -2
    set max = 2
  endif
  set delta = `echo "scale=2; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba${atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk 'NR>1 {print $1, log($6)/log(10)}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif


####################
# Nusselt number
####################
if ($nusselt == 1) then
  set adaptive_scale = 1
  echo "Nusselt number time series"

  set title = ""
  set xlabel = "Time"
  set zlabel = "Nusselt_number"

  set pdfoutput = nusselt_ts.pdf
  set jpgoutput = nusselt_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk 'NR > 1 {if(min==""){min=$2}; if($2< min) {min=$2}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk 'NR > 1 {if(min==""){min=max=$2}; if($2>max) {max=$2}} END {printf "%12.5f", max}' $timeseries`
  else  
    set min = 0
    set max = 20
  endif
  set delta = `echo "scale=1; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba{$atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk 'NR>1 {print $1, $2}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

########################
# Viscous dissipation
########################
if ($viscdiss == 1) then
  set adaptive_scale = 1
  echo "Viscous dissipation time series"

  set title = ""
  set xlabel = "Time"
  set zlabel = "log@-10@-(Viscous_dissipation)"

  set pdfoutput = viscdiss_ts.pdf
  set jpgoutput = viscdiss_ts.jpg

  set min = `awk 'NR > 1 {if(min==""){min=log($7)/log(10)}; if(log($7)/log(10)< min) {min=log($7)/log(10)}} END {printf "%12.5f", min}' $timeseries`
  set max = `awk 'NR > 1 {if(min==""){min=max=log($7)/log(10)}; if(log($7)/log(10)>max) {max=log($7)/log(10)}} END {printf "%12.5f", max}' $timeseries`
  set delta = `echo "scale=1; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba${atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk 'NR>1 {print $1, log($7)/log(10)}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

########################
# Composition
########################
if ($composition == 1) then
  set adaptive_scale = 1
  echo "Composition time series"

  set title = ""
  set xlabel = "Time"
  set zlabel = "Composition"

  set pdfoutput = composition_ts.pdf
  set jpgoutput = composition_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk 'NR > 1 {if(min==""){min=$8}; if($8< min) {min=$8}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk 'NR > 1 {if(min==""){min=max=$8}; if($8>max) {max=$8}} END {printf "%12.5f", max}' $timeseries`
    if ($min == $max) then
      set min = -1
      set max = 1
    endif
  else
    set min = 0.2 
    set max = 0.8
  endif
  set delta = `echo "scale=4; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba${atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk 'NR>1 {print $1, $8}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

########################
# Internal heating
########################
if ($heats == 1) then
  set adaptive_scale = 1
  echo "Internal heating time series"

  set title = ""
  set xlabel = "Time"
  set zlabel = "RaQ"

  set pdfoutput = heating_ts.pdf
  set jpgoutput = heating_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk 'NR > 1 {if(min==""){min=$17}; if($17< min) {min=$17}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk 'NR > 1 {if(min==""){min=max=$17}; if($17>max) {max=$17}} END {printf "%12.5f", max}' $timeseries`
    if ($min == $max) then
      set min = -1
      set max = 1
    endif
  else
    set min = 0.2 
    set max = 0.8
  endif
  set delta = `echo "scale=4; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba${atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk 'NR>1 {print $1, $17}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

########################
# Mixed fraction
########################
if ($mixfrac == 1) then
  set adaptive_scale = 0
  echo "Mixed fraction time series"

  set title = ""
  set xlabel = "Time"
  set zlabel = "Mixed_fraction"

  set pdfoutput = mixfrac_ts.pdf
  set jpgoutput = mixfrac_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk 'NR > 1 {if(min==""){min=$16}; if($16< min) {min=$16}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk 'NR > 1 {if(min==""){min=max=$16}; if($16>max) {max=$16}} END {printf "%12.5f", max}' $timeseries`
    if ($min == $max) then
      set min = -1
      set max = 1
    endif
  else
    set min = 0
    set max = 1
  endif
  set delta = `echo "scale=4; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba${atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk 'NR>1 {print $1, $16}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

########################
# Depletion
########################
if ($depletion == 1) then
  set adaptive_scale = 0
  echo "Depletion time series"

  set title = ""
  set xlabel = "Time"
  set zlabel = "Depletion"

  set pdfoutput = depletion_ts.pdf
  set jpgoutput = depletion_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk 'NR > 1 {if(min==""){min=$18}; if($18< min) {min=$18}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk 'NR > 1 {if(min==""){min=max=$18}; if($18>max) {max=$18}} END {printf "%12.5f", max}' $timeseries`
    if ($min == $max) then
      set min = -1
      set max = 1
    endif
  else
    set min = 0
    set max = 1
  endif
  set delta = `echo "scale=4; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba${atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk 'NR>1 {print $1, $18}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

########################
# Stagnant lid
########################
if ($staglid == 1) then
  set adaptive_scale = 0
  echo "Stagnant-lid time series"

  set title = ""
  set xlabel = "Time"
  set zlabel = "Stagnant_lid"

  set pdfoutput = staglid_ts.pdf
  set jpgoutput = staglid_ts.jpg

  if ($adaptive_scale == 1) then
    set min = `awk 'NR > 1 {if(min==""){min=$19}; if($19< min) {min=$19}} END {printf "%12.5f", min}' $timeseries`
    set max = `awk 'NR > 1 {if(min==""){min=max=$19}; if($19>max) {max=$19}} END {printf "%12.5f", max}' $timeseries`
    if ($min == $max) then
      set min = -1
      set max = 1
    endif
  else
    set min = 0
    set max = 1
  endif
  set delta = `echo "scale=4; ($max - $min)/5" | bc -l`

  psbasemap -X$X0 -Y$Y0 -R$time0/$time1/$min/$max/ -JX$dimx/$dimz -Ba${atime}:${xlabel}:/a{$delta}:${zlabel}:WSen:.${title}: -K >! $output
  awk 'NR>1 {print $1, $19}' $timeseries | \
  psxy -R -JX -W5/0/0/0 -O >> $output

  ps2raster $output -Au -P -Tf -F$pdfoutput
  if ($makejpg == 1) then
    ps2raster $output -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  endif

endif

##############################
# Moving and cleaning
##############################

if ($makejpg == 1) then
  mv -f *.jpg $dirplots
endif

if ($makejpg == 1 && $deletepdf == 1) then
  rm -f *.pdf
else 
  mv -f *.pdf $dirplots
endif

rm -f *.ps 
